
clear all;
close all

%%%%%%%%%% Don't change the above code %%%%%%%%%%%%
global theta1 dbar  L gL alpha rbest rho sigma1 ns n betas d delta1 betad

t1 = 0.2;
t2 = 0.8;
nt = 100;
staux = t1:(t2-t1)/(nt-1):t2;
staux = 0.648;
nt=1;
for itaux=1:nt

    taux=staux(itaux);
    taf = 0;

    n     =  2; % number of countries
    ns      =  1;
    betad   = 1/(1+0.1);     % discount factor
    theta1  = 4;
    sigma1 = 2;
    gL      = 0.02;
    dbar   = 1.4;
    delta1 = 1-1/(1+gL);


    L = ones(n,1);
    betas = ones(ns,1)/ns;

    rho   = betad*(1-delta1);
    rbest = delta1/(theta1*(1-betad*(1-delta1))+delta1);

    ncountries=n;
    nsectors = ns;


    alpha = [1;0.9];

    % private eqm
    r_private = rbest;
    n  = ncountries;
    ns = nsectors;
    ilog_rL = 1; % use exp(rL)/(1+exp(rL))

    dij = dbar;
    d=dij*ones(n); %dij is (1+iceberg) cost producing in j selling to i
    d=d-(dij-1)*eye(n);


    %% solve for the steady state of gov assuming dG2/dT=0!!!
    % This is not the true steady state, it's only for us to get initial guess
    % --- steady state ---


    options = optimoptions('fsolve', 'TolFun', 1e-12);
    options.Display = 'iter';

    xguess(1:2) = 1.1;



    [xval,fval_sol,eflag_ss] = fsolve(@(xval) solve_1s_2c_private_notForRamsey(xval,taux,taf),xguess,options); %iter


    w(1:n,1) = xval(1:n); % row vector
    taf = 0;

    r1 = log(rbest);
    rL(1,1) = exp(r1)*L(1);%rbest*L(1)
    rL(2,1) = rbest*L(2);

    T = alpha.*rL/delta1;

    sextax = zeros(n);
    for i=1:n-1
        sextax(i+1) = taux;
    end
    staf = zeros(n);

    for i=1:n-1
        staf(i+1) = taf;
    end

    pi = zeros(n,n); % (1,2,3): country 2 export to ecountry 1 in sector 3, NOTE country order is opposite to our note
    %d, matrix of n times n, dij country j ship to country i

    xnimatrix1= T'.*((w'.*(1+sextax).*(1+staf).*d).^(-theta1));
    xn1 =sum(xnimatrix1,2).*ones(n,n);

    pimatrix1 = xnimatrix1./xn1;
    pi = pimatrix1;
    xnimatrix= xnimatrix1;

    betas =1;
    x  =  (1+theta1)/theta1*w.*(L-sum(rL,2));

    for m=2:n
        x(1)=x(1)+betas(1)*taux/(1+taux)*pi(m,1)*x(m);
    end
    tempt=0;

    for m=2:n
        tempt=tempt+betas(1)*taf/(1+taf)*pi(1,m);
    end

    x(1)=x(1)/(1-tempt);

    P=(sum(xnimatrix1,2)).^(-betas./theta1);

    fval(1) = P(1)-1;
    fval(2) = pi(1,2)*x(1)-pi(2,1)*x(2);

    Tss = T;
    xss = x;
    Pss =P;
    pi_ss=pi;



    xval_ss= xval;
    wss = xval(1:n);
    taf_ss = 0;
    rss = rbest;

    sx(:,itaux) = x;
    sT(:,itaux) = T;
end

